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ABSTRACT 

Physical collisions between stars occur frequently in dense star clusters, either via 
close encounters between two single stars, or during strong dynamical interactions 
involving binary stars. Here we study stellar collisions that occur during binary-single 
and binary-binary interactions, by performing numerical scattering experiments. Our 
results include cross sections, branching ratios, and sample distributions of parameters 
for various outcomes. For interactions of hard binaries containing main-sequence stars, 
we find that the normalized cross section for at least one collision to occur (between 
any two of the four stars involved) is essentially unity, and that the probability of 
collisions involving more than two stars is significant. Hydrodynamic calculations have 
shown that the effective radius of a collision product can be 2-30 times larger than 
the normal main-sequence radius for a star of the same total mass. We study the 
effect of this expansion, and find that it increases the probability of further collisions 
considerably. We discuss these results in the context of recent observations of blue 
stragglers in globular clusters with masses exceeding twice the main-sequence turnoff 
mass. We also present Fewbody, a new, freely available numerical toolkit for simulating 
small- iV gravitational dynamics that is particularly suited to performing scattering 
experiments. 

Key words: stellar dynamics - methods: A'^-body simulations - methods: numerical 
- binaries: close - blue stragglers - globular clusters: general. 



1 INTRODUCTION 

Close encounters and direct physical collisions between stars occur frequently in globular clusters. For a star in a dense cluster 
core, the typical colli sion time can be c omparable to the cluster lifetime, implying that essentially all stars could have been 
affected by collisions fc Daviri976ll . E ven in moderate ly dense clusters, col l isions can hap pen frequently during resonant 

inter a ctions involving primordial binaries jHut fc Verbunt| [l983: L conardlllQsl ISieurdsson fc Phinncy 1991; IPavies fc Ben3 
Il99,4 Is'jgurdsson fc Phinnevlll99!Tl : iBacon et alJll99fi^ . In open clusters with significant binary fractions (~ 10% or more), 
mergers may occur mo re often through binary-b inary interactions than through single-single collisions and binary-single 
interactions combined jbeonard fc FahlmarJ Il99lll . Collisions involving more than two stars can be quite common during 
binary-single and binary-binary interactions, since the product of a first collision between two stars expands adiabatically 
following shock heating, and therefore has a larger cross section for subsequent collisions with the remaining star(s). 

Collisions and binary interactions strongly affect the dynamical evolution of globular clusters. The formation of more 
massive objects through mergers tends to accelerate core collapse, shortening cluster lifetimes. On the other hand, mass 
loss from evolving collision products can indirectly heat the cluster core, thereby postponing core collapse. The realization 
during the 1990s that primordial binaries are present in globular clusters in dynamically significant numbers has completely 
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chan ged our theoretical perspective on these systems iGoodman fc Hulll989l : ISieurdsson fc PhinnevI FlQQSl : llvanova et alJ 
[200^. Most importantly, dynamical interactions of hard primordial binaries with sin gle stars and other binaries are thought 



IZUU4II . Most importantly, dynamical interactions ot nard primordial binaries with sin gle stars and other pmaries are thouRht 
to be the primary mechanism for supporting globular clusters against core collapse (M£MillaiL_eLjyJ[l993, [iMl^ 
IQmllHut et alllflflilHeggie Aarsethll 99d:lMcMil]a,n HuiIi994Ir^o et al.l2nnil:lFregea,u et al .1200.4 lOiersz Snurzen- 



I 



20Qj ) . Observational evidence for the existence of primordial binaries in globular clusters is now well established ijHut et al. 



I992I: ICote et al.lfl994l : iRubenstein fc Bailvnlll997il . Recent Hubble Space Telescope {HST) observations have provided direct 
constraints on the primordial binary fractions in many clusters. For example, the observation of a broadened main sequence in 
NGC 6752, based on HST-'WFP C2 images, suggests that t he binary fraction is pr obably in the range 15%-40% in the inner 
core jRubenstein fc Bailvnlll997ll . Using a similar method. iBellazzini et al.l (I2OO2I1 find that the binary fraction in the inner 
region of NGC 288 is probably between 10% and 20%, and less than 10% in the outer region. Observations of eclip sing binaries 
and BY Draconis stars in 47 Tuc yield an estimate of ~ 13% for the core binary fraction llAlbrow et ahlEoOlh . although a 
recent reinterpr etation of the observ ations in combination with new theoretica l results suggests that this number might be 
closer to ~ 5% Jivanova et alJl2004^ . Using HST-WFPC2. iBolton et al.l Jl999l) derive an upper limit of only ~ 3% on the 
binary fraction in the core of NGC 6397. 

In this paper, we focus on interactions involving main-sequence (hereafter MS) stars, and the production of blue stragglers 
(hereafter BSs). BS stars appear along an extension of the MS blue- ward of the turnoff point in the colour-magnitude 
diagram (CMD) of a star cluster. All observations suggest t hat they are massive MS stars formed through mergers of two (or 
more) lower-mass stars. For example. ICilliland et al.l Jl998l) have demonstrated that the masses estimated from the pulsation 
frequencies of four oscillating BSs in 47 Tuc are consistent with their positions in the CMD. Indirect measurements of BS 
masses yield values of up to four tim es the MS turnoff mass, although the uncertainties are significant ([Bond fc Perrvlll97ll : 
IStrom et al]ll97ll: iMilone et a'illl992l ) . More recent spectroscop ic measurements yield much more precise masses, with one BS 
in 47 Tuc about twic e the MS turnoff mass JShara et al.lll997l) . and two in NGC 6397 more than twice the MS turnoff mass 
JSepinskv et al.ll2000ll . 

Mergers of MS stars c an occur in at l east two different ways: via physical collisions, or througjh the coalescence of two stars 

1 ll98fll : lLiviolll993l: Istrvkeilll993 : iBailvn fc Pinsonneaultlll99a . 



in a close binary system l|Leonardlll98fll : lLiviolll993l : IStrvkeiill99.1 : IBailvn fc Pinsonneaiiltlll99Hfl . Direct evi dence for binar y 
progenitors has been found in the form of contact (W UMa t ype) binaries among B Ss in many glo bular clusters l|^icin ski200(]ll . 
including low-densit y globular clusters such as NGC 288 (iBeUazzini et a l.' 2002), NGC 5466 iUateo et al.lll990ll. and M71 
jYan fc Mateclll994l) . as well as in many open clusters fe.g.. ljaln^^dTIl995. : .Kaluznv fc Rucin^klE^g^nMUone fc LathamI 
Il994^ . At the same time, strong indication for a coUisional origin comes fror n detections by HST of large numbers of bright 
BSs concentrated in the c ores of some of the densest clusters, such as M15 Jde Marchi fc Parescelll994l: IVannv et alJll994al : 
ICuhathakurta et aljll996l) . M30 JVannv et alJll994bl: Ic^hathakurta et al.lll99a^ . NGC 6397 ijBurgarella et al.lll994E . NGC 
6624 JSosin fc Kin J M95rl . and M80 ^Fe^^^^^l.lll99S^ . High-resolution HST images reveal that the central density profiles 
in many of these clusters steadily increase down to a radius of ~ 0.1 pc, with no signs of flattening. Direct stellar collisions 
should be extremely frequent in such high density environments. 

Evidence for the greater importance of binary interactions over direct collisions of single stars for producing BSs in 
some globular clusters can be found in a lack of correlation between BS specific frequency and cluster central collision rate 
ide Angeh fc Piott dboosi: iFerraro et alJbooj)^ More direct evidenc e comes from the BS S1082 in the open cluster M67, which 
is part of a wide hierarchical triple svstem iSandauist et aljl2003l) . The most natural formation mechanism is via a binary- 
binary interaction. There is further evidence in the radial distributions of BSs in clusters. HST observations, in combination 
with ground-based studies, have revealed that the radial distributions of BSs in the clusters M3 and 47 Tuc are bimodal - 
peaked in the core, decreasing at intermediate radn, and rising again at larger radn JFerra.ro et al J B F^. The 

most plausible explan ation is that the BSs a t larger radii were formed through binary interactions in the cluster core and 
ejected to larger radii dSigurdsson et aDll994^ . 

In this paper, we perform numerical scattering experiments to study stellar collisions that occur during binary interactions. 
One approach for attacking the problem is to perform a full globular cluster simulation, taking into account every relevant 
physical process, including stellar dynamics, stellar evolution, and hydrodynamics. This approach is enticing in its depth, but 
would certainly yield results with a complicated dependence on the input parameters and physics that would be difficult to 
disentangle. A simpler approach is to study in detail the scattering interactions that occur between binaries and single stars or 
other binaries. This approach isolates the relevant physics and produces results that are easier to interpret. Furthermore, the 
cross sections tabulated will be useful for future analytical and numerical calculations of cluster evolu tion and interaction rates. 
For a discussion of the interplay between globular cluster dynamics and stellar collisions, see, e.g.. iHurlev et all teoOlh . For 
dense globular cluste r cores, merger rates via binary stellar evolution can be significantly enhanced by dynamical interactions 
Jivanova et alJl2004^ . 

Our paper is organized as follows. In § 2 we summarize previous theoretical work on stellar collisions in binary interactions. 
In § 3 we describe our numerical method, and introduce the two numerical codes used. In § 4 we test the validity of our 
numerical method by comparing with previous results. In § 5 we present a systematic study of the dependence of the collision 
cross section in binary-single and binary-binary interactions on several physically relevant parameters. In § 6 we consider 
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binaries with parameters characteristic of those found in globular clusters, and study the properties of the resulting binaries 
and triples containing collision products. Finally, in § 7 we summarize and conclude. 



2 PREVIOUS WORK 

The re now exists a very l arge body of numerical work on bi nary-single and , to a l esser extent, binary-binary interactions (see, 
e.g., iHeggie fc Hu^booi for an overview and references) . iHut fc Bahcalll lll983ll perfo rmed one of the most extensive early 
studies of binary-single star scattering for the equal-mass, point-particle case. iMikkoS Hiiii) performed the first systematic 
studies o f binary-binary interactions in the point-particle limit, first for the case of equal energy binaries, and later for unequal 
energies (lMikkolalll984al). He also studied t he energy generated in binary-binary interactions in the context of the evolution 
of globular clusters jMikkolallTgSSbl Il984bh . Most numerical scattering experiments have been performed in the point-mass 
limit, neglecting altogether the effects of the finite size of stars. However, as we summarize below, there are a number of studies 
that apply approximate prescriptions for dissipative effects and collisions post facto to numerical integrations performed in the 
point-mass limit but in which pairwise closest approach distances were recorded. There is also one study in which collisions 
are treated in situ in a simplified manner, and several that perform full smoothed particle hydrodynamics (SPH) simulations 
of multiple - star in teractions (see below). 

iHoffej ijlflSlj ^ was the first to study distances of closest approach between stars in both binary-single and binary-binary 
interactions. He found that roughly 40% of bin ary-binary encounters in a typical globular cluster core will lead to a physical 
collision between two stars. iKrolik et all il984h considered the evolution of a compact binary in a globular cluster core subject 
to perturbations by single field sta r s, and found that an induced merger or collision between two stars in a binary-single 
interaction is likely. iHut fc Inaeakil l|l98.4 l applied a single-parameter, "fully inelastic sphere" collision model after the fact 
to a l arge number of b inary-single interactions for which distances of closest approach were recorded, and calculated merger 
rates . iMcMillanI (Il98dl applied simple prescriptions for the dissipative effects of gravitational radiation, tidal interactions, and 
physical contact between stars after the fact to a large number of binary-single interactions involving tight binaries. He found 
that dissipative effects reduce binary heating efficiency in cluster cores by roughly an order of magnitude over that obtained in 
the point-mass limit. He also found that the most likely outcome of a bin ary-single interaction involving a tight binary is the 
coalescence of at least two of the stars. This work was carried further bv lPortegies Zwart et all l|l997al lb[l. who included the 
effects of binary stellar evolution on binary-single interactions. They found that about 20 per cent of encounters between a 
primordial b inary and a clus ter star result in collisions, while almost 60 per cent of encounters with tidal-capture binaries lead 
to collisions. iLeonarcJ (Il989l) performed a small nu mber of binary-binary interactions and recorded close approach distances 
and calculated ejection speeds of collision products. iLeonard fc FahlmanI l)l99lll performed the first set of binary-single and 
binary-binary interactions in which stars were allowed to merge during the numerical integration. They studied the rate 
of productio n of BSs in clusters, and performed the first, simplified "population synthesis" study of BSs in clusters. iHillj 
il99ll Il992f) considered stars with a range of masses exchanging into binaries, and found the distance of closest approach 
to be roughly a constant fr action of binary semimajor axis independent of intruder mass, over a wide range of mass ratios. 
Iciearv fc MonaghanI l)l99nl) performed full SPH si mulations of binary-single inter actio ns and showed direct l y the importance 
of taking into account the non-zero size of stars. lOoodman fc HernauistI l)l99lh and IPavies et al.l l|l99.'j . Il994h performed 
sets of binary-single and binary-binary interactions with tight binaries in t he point-mass limit , and selected a handful to 
run with a full SPH code. They found that multiple mergers are common. iBacon et al.l il996tl performed a large number 
of binary-binary interactions and presented a survey of close approach cross sections for several sets of physically relevant 
binary parameters. They also calculated outcome frequencies, studied the properties of the inter action products, and used 
their results in analytical calculations of interaction rates in globular cluster cores. More recently, lOiersz fc SpurzemI l|200,'j) 
have incorporated into their Monte-Carlo globular cluster evolution code Aarseth's NBODY iox performing direct integrations 
of binary interactions. By storing the results of binary interactions that occur during cluster evolution, they have calculated 
close approach cross sections, and a few differential cross sections. 



3 NUMERICAL METHOD 

The scattering experiments presented in this paper were performed primarily using Fewbody, a new numerical toolkit designed 
for simulating small- A'' gravitational dynamics, which we describ e below. In some cases we also use the scattering facilities 
of the Starlab software environment l|Portegies Zwart et alJl200l|) . Starlab was used mainly to compare with Fewbody, but in 
cases where Starlab data were compiled before Fewbody was written, Starlab results were used. In particular, all calculations 
in §|S|were performed with Starlab. 
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3.1 Setup 

We label the two objects in a scattering experiment and 1. In the case of binary-single scattering, is the binary and 1 is 
the single star. In the case of binary-binary scattering, and 1 are each binaries. We use the same system of labelling for each 
binary, so the members of binary i are labelled iQ and il. There are several parameters required to uniquely specify a binary- 
single or binary-binary scattering experiment. To describe the initial hyperbolic (or parabolic) orbit between objects and 1, 
one needs to specify the relative velocity at infinity Voo , impact parameter b, and masses mo and mi . To describe the internal 
properties of each object, one needs to specify the semimajor axes a;, eccentricities d, individual masses rriij, and stellar 
radii Rij. There are also several phase and orientation angles required for each binary: the orientation of the binary angular 
momentum vector relative to the angular momentum vector describing the orbit between and 1, given by the polar angle 9 
and the azimuthal angle 0; the angle ui between the binary Runge-Lenz vector and some fiducial vector perpendicular to the 
binary angular momentum vector (e.g., the cross product of the binary angular momentum and the 0-1 angular momentum); 
and ?7, the mean anomaly of the binary. For all the scattering experiments presented in this paper, these phase and orientation 
angles are chosen randomly, so that the cross sections calculated represent averages over these quantities. In detail, these 
angles are given hy 9 — cos~^(2X — 1), </> = 2nX, uj = 2-kX, and rj = 2-kX, where X is a uniform deviat e in the ran ge [0, 1). In 
addition, unless otherwise noted, the eccentricity of each binary is chosen from a thermal distribution l^,Ieanjll919^ truncated 
at large e such that there is no contact binary. In each scattering experiment, numerical integration is started at the point at 
which the tidal perturbation (-Ftid/^rei) on a binary in the system reaches 5 (see 13. 3. 411 . 

It is customary to specify the relative velocity at infinity in terms of the critical velocity, Vc, defined such that the total 
energy of the binary-single or binary-binary system is zero. For Voa > Vc the total energy of the system is positive, and full 
ionization is possible. That is, a possible outcome of the scattering experiment is that each star leaves the system unbound 
from any other with positive velocity at infinity. For Voc < Vc, the total energy of system is negative, and the encounters are 
likely to be resonant, with all stars involved remaining in a small volume for many dynamical time-scales. Defining the total 
mass M — mo + mi and reduced mass /i = momi/M, the critical velocity is 

1/2 

(1) 



(2) 





/ moomol \ 




\ ao J 



for the binary-single case, and 

G / moompi ^ miomn 
fi \ ao ai 



for the binary-binary case. The cross section for outcome X is obtained by performing many scattering experiments out to a 
maximum impact parameter 6max and calculating 

where Nx is the number of experiments that have outcome X, and A'' is the total number of scattering experiments performed. 
In all cases the maximum impact parameter was chosen large enough to ensure that the full region of interest was sampled. 
In other words, for 6 > &max, all interactions are fiy-by's in which each binary is only weakly tidally perturbed during the 
interaction. For calculations performed with Fewbody, &max was chosen to correspond to a pericentre distance of rp = 5(oo-|-ai) 
in the binary-binary case, and Vp = 5a in the binary-single case. For this value of pericentre distance, the binary eccentricity 
induced in the fiy-by is quite small (5e ^ 1 for initially circular binaries, and 5e/e <^ 1 for non-circular binaries; see 
iRasio fc Heggielll99 st iHeggie fc Rasiolll996l) . For calculations performed with Starla b, &max was chosen auto matically by using 
successively larger impact parameter annuli until no relevant outcomes were found jMcMillan fc HuJl996t) . The uncertainty 
in the cross section is calculated assuming Poisson counting statistics, so that 



Acta- = Ti-femax ^ ■ (4) 

In principle, it is nec essary to include scattering experiments that result in unresolved outcomes in this uncertainty (see, e.g., 
iHut fc Bahcalil983l) . However, in practice we find that the number of unresolved outcomes is small, and does not significantly 
contribute to Acta. 



3.2 Possible Outcomes 

The possible outcomes of binary-single and binary-binary scattering interactions are listed in Tables Q and |21 ordered by 
the number of collisions, ricoii. Stars are represented as filled circles, brackets enclose two objects that are bound to each 
other in a binary, and colons represent physical collision products. In Table|2|we also list the abbreviations used in the paper 
to refer to certain outcomes. When there are no collisions (as is the case in the point-mass limit), the number of possible 
outcomes is small, as shown in the ncoii = rows in each table. However, when one considers stars with non-zero radius 
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and allows for the possibility of collisions and subsequent mergers, the total number of outcomes becomes large. Assuming 
indistinguishable stars, there are six possible outcomes for the binary-single case, and 16 for the binary-binary case. These 
numbers are evidently increased for distinguishable stars. The software used in this paper distinguishes among all possible 
outcomes. 



3.3 Fewbody 

Fewbody is a new numerical toolkit for simulating small- gravitational dynamics. It is a general A-body dynamics code, 
although it was written for the purpose of performing scattering experiments, and therefore has several features that make it 
well-suited for this purpose. It can be described succinctly in terms of its key elements. 



3.3.1 Adaptive Integration and Regularization 

At its core, Fewbody uses the 8th-order Runge-Kutta Prince-Dormand integration method with 9th-order error estimate 
and adaptive timestep to advance the A-body system forward in time. It integrates the usual formulation of the Af-body 
equations in configuration s pace, but allows for the option of global pairwise Kustaanheimo-Stiefel (K-S) regularization 
l)Heggi(Jll97j : iMikkolalHoSsh . Global regularization is a coordinate transformation that removes all singularities from the 
A^-body equations, making the integration of close approaches, and even collision orbits, much more accurate. It is well-suited 
for small- A' dynamics, since it requires the integration of ~ A^ separations instead of A^ positions, and becomes prohibitively 
computationally expensive for A^ > 10. Although it should in principle make numerical integration more accurate, it was found 
that the adaptive timestep algorithm alone performed as well as global regularization, in terms of the computational time 
required for a specified level of energy and angular momentum accuracy. The use of regularization requires extra effort to 
detect physical collisions, since, with regularization, pericentre is not necessarily resolved by the integrator. For the sake of 
si mplicity, we ha ve chosen not to implement the appropriate technique for detecting collisions with regularization (see § 9.8 
of lAarsethir2003t) . Furthermore, physical collisions naturally soften the singularities in the non-regularized A'-body equations 
by making them physically inaccessible. Regularization was therefore only used to test calculations made in the point-mass 
limit. For all other calculations the non-regularized integration routine was used. 



3.3.2 Classification 

Fewbody uses a binary tree algorithm to handle several aspects related to performing scattering experiments. Most importantly, 
it uses a binary tree algorithm to classify the A'-body system into a set of independently bound hierarchies. For example, if 
the outcome of a scattering experiment between two hierarchical triples is a hierarchical triple composed of binaries, Fewbody 
will classify it accordingly. Fewbody creates the set of binary trees iteratively, according to the following simple rules. First, 
as shown in Figure any existing set of trees is flattened so that each star in the A'-body system represents the top-level 
node of a one- node tree. Next, as shown in Figure |5| the two top-level nodes that are bound to each other with the smallest 
semimajor axis are replaced by a parent node containing all dynamical information about the centre of mass, as well as all 
information about the binary's orbit, including phase. The previous step is repeated, as shown in Figure^! until no top-level 
nodes are found to be bound to each other. This algorithm is clearly general in A. The resulting set of binary trees is a unique 
classification of the configuration of the A'-body system. As described below, the classification is used for determining when 
an interaction is complete. The binary tree algorithm is also used (with a slightly different set of rules for creating the trees) 
to make the numerical integration more efficient, as also described below. 



3.3.3 Stability 

Fewbody assesses the dynamical stability of gravitationally bound hierarchies in an approximate way using the classification just 
described, and a simple analytical test. There currently exists only one reasona bly accurate criterion f or the dynamical stability 
of an A > 2 gravitational system, the approximate analytical criterion of iMardline fc AarsethI l|200lh for the dynamical 
stability of hierarchical triples. Fewbody assesses the stability of each binary tree by applying this criterion at each level in the 
tree. For example, for a hierarchical quadruple system (which consists of a star in orbit around a hierarchical triple - shown 
as "[[[• •] •] •]" in Table |5J, it first applies the triple stability criterion to the inner triple, then applies it to the "outer" 
triple, treating the innermost binary as a single object. For the case of a hierarchical quadru ple composed of two b inarie s 
("[[• •] [• •]]" in Table|2J, Fewbody uses the additional correction factor presented in § 4.2 of lMardling fc AarsethI (|20fllh . 
The stability of a hierarchical system as determined by this method is only approximate, but, in our experience, seems to 
work reasonably well. For the particular case of binary-binary scattering, hierarchical triples which appear to be stable are 
classified as unstable less than roughly one percent of the time. 

It should be noted that the stability assessed here is dynamical rather than secular, so, e.g., any resonances that would 



© 2004 RAS, MNRAS 000, 000-000 



6 J. M. Fregeau, P. Cheung, S. F. Portegies Zwart, & F. A. Rasio 



destroy an otherwise stable hierarchical system are ignored. Such resonances are likely to be important in the more genera l 
context of the dynamics of globular clusters and their constituent populations fe.g.. lMiller fc Hamiltorl2002l : librd et al.l200d) . 
but are beyond the scope of the present paper. It should also be noted that our use of binary trees prevents us from recognizing 
stable three-body systems which are not hierarchical, such as the stable "figure eight" orbit for three stars of comparable mass 
ijchenciner fc MontgomervlboOfi : iMontgomervl l200lh . Howeve r, the fraction of strong binary-binary scattering encounters 
resulting in this configuration is likely to be exceedingly small l)Heggiell200(il ). 



3.3.4 Hierarchy Isolation 

Fewbody also uses a binary tree algorithm to speed up numerical integration by isolating from the integrator certain tight 
binaries and hierarchies that are only very weakly perturbed, yet dominate the calculation by driving down the integration 
time-scale. It does this by integrating only the top-level nodes (centres of mass) of a set of binary trees created using the 
algorithm described in S 18.3.21 but subject to a slightly different rule-set. Two top-level nodes can only be replaced by their 
parent node if: (1) the binary tree represented by the parent node is stable, (2) the tidal perturbation on the outer binary of 
the tree (the two nodes below the top-level) at apocentre due to all other top-level nodes in the system is less than a specified 
fraction, S, of the minimum force between them (-Ftid/^roi < S at apocentre), and (3) the evolution of the binary tree can be 
treated analytically. The relative force at apocentre is calculated simply as 

Gmomi , , 

[„(l + e)p' 

where mo and mi are the masses of the members of the outer binary, a is the semimajor axis, and e is the eccentricity. The 
tidal force at apocentre is calculated simply as 

2G(mo -f mi)mi , 
Ftid = 2_^ ^ -3 a{l + e), (6) 

i 

where the sum is taken over all other top-level nodes in the system, m^ is the mass of the other top-level node, and is the 
distance to the other top-level node. Note that this sum represents the upper limit of the tidal force since it does not take 
into account relative inclination between the binary and the other top-level nodes. 

A binary (or hierarchy) that is isolated from the integrator in this way is treated numerically again when its relative 
tidal perturbation exceeds 5. This is done by resuming the integration from the previous step (when the hierarchy's tidal 
perturbation was less than 5) with the parent node replaced by its child nodes, and orbital phase advanced to the current time. 
In practice, this algorithm isolates from the integrator mainly weakly-perturbed binaries, and a few extremely hierarchical 
triples in which the tidal perturbation on the inner binary due to the outer member is very small. For binary-single scattering, 
hierarchy isolation can speed up the integrations by up to an order of magnitude on average. For binary-binary scattering, 
especially when the two binaries have very disparate semimajor axes, and hence orbital time-scales, this algorithm can speed 
up the integrations by a few orders of magnitude on average. The quantity 5 plays the role of an integration tolerance 
parameter. Larger values of 5 allow hierarchies to be treated analytically more frequently, yielding faster calculations but 
sacrificing energy accuracy. Smaller values of 5 yield better energy conservation at the expense of computational speed. 



3.3.5 Calculation Termination 

Fewbody uses the classification and stability assessment techniques outlined above, in combination with a few simple rules to 
automatically terminate the integration of scattering encounters when they are complete - in other words, when the separately 
bound hierarchies comprising the system will no longer interact with each other or evolve internally. Integration is terminated 
when: (1) each pair of top-level nodes has positive relative velocity, (2) the tidal perturbation (-Ftid/^rei) on the outer binary 
(the two nodes below the top-level) of each tree due to the other top-level nodes is smaller than S, (3) each tree is dynamically 
stable (as defined in § I3.3.3II . and (4) the A''-body system composed of the top-level nodes has positive energy. The last 
condition is required because it is possible for the members of an A*'-body system to be separately unbound and receding from 
each other, yet for the system as a whole to be bound. Here 5 again plays the role of an accuracy parameter, with smaller S 
yielding more accurate outcome classifications. Since the A''-body problem is chaotic, with initially neighbouring trajectories 
in phase space diverging exponentially, the value of S should play only a minor role in the statistical accuracy of classifications 
of outcomes. 



3.3.6 Physical Collisions 

Fewbody performs collisions between stars in the "sticky star" approximation. In this approximation, stars are treated as 
rigid spheres with radii equal to their stellar radii. When two stars touch, they are merged with no mass lost, and with 
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linear momentu m conserved. (Tidal effects, which may significantly increase the collision rate for close encounters (see, e.g., 
iMcMillanlEgiil . are beyond the scope of this method, but may be approximated by larger initial effective stellar radii.) The 
radius of the merger product is set to 

^merger = fcxp{Rl + ^2) , (7) 

where Ri and R2 are the radii of the merging stars, and /exp is an expansion factor. To determine a reasonable value for /cxp, 
one must consider the relevant time-scales involved. The characteristic time-scale of a typical binary scattering encounter in 
a globular cluster core is between ~ 10 yr for a fly-by and < 10* yr for a resonant encounter, while the thermal time-scale of a 
~ i Mq MS star is ~ 10^ yr. Therefore, it is invalid to treat merger products as rejuvenated ("reborn") MS stars (/cxp — 1) 
during scattering encounters. The hydrodynamical time-scale is ~ 1 hr, so it is more accurate to treat merger products as 
hydrodynamically settled. SPH simulations show that /cxp should be in the range 2-30, depending on the relative orientations 
of the two stars before collision l|Lombardi et aljbofl.'j) . These simulations also show that the amount of mass lost in the types 
of collisions characteristic of globular clusters is typically of order 1%, so our assumption of zero mass loss is a reasonable first 
approximation . 

Collision products are likely to have significant rotation and be non-spherical. Furthermore, it is not clear that the value 
of the expansion parameter for the merger of two pristine MS stars should be the same as that for mergers involving collision 
products. Thus /cxp should be considered an effective quantity, averaged over many collisions. A more realistic approach that 
adopts several separate parameters is in principle possible, but beyond the scope of the current paper. 



3.3.7 General Availability 

Fewbody is freely available for download on the web^, licensed under the GNU General Public License (GPL). It contains a 
collection of command line utilities that can be used to perform individual scattering and A^-body interactions, but is more 
generally a library of functions that can be used from within other codes. Its facilities make it aptly suited for performing 
scattering interactions from within larger numerical codes that, e.g., calculate cross sections, or evolve globular clusters via 
Monte-Carlo techniques. 

Available along with Fewbody, there is an OpenGL-based visualization tool called GLStarView that can be used to view 
A'^-body interactions as they are being calculated by Fewbody, in an immersive, 3-D environment. GLStarView has proven to 
be a valuable aid in developing our understanding and physical intuition of binary interactions. 



3.4 Starlab 

Starlab is a collect ion of modular software too ls designed to simulate the evolution of dense stellar systems and analyse the 
resulting data (see IPorteeies Zwart et al]l200ll . for a detailed description). It is freely available on the web^. It consists of a 
library of programs for performing stellar dynamics, stellar evolution, and hydrodynamics, together with a set of programs 
acting as bridges between them. They may be combined to study all aspects of the evolution of A-body systems. For this 
paper, we use the three-body scattering facility scatter3 and the general A-body scattering facility scatter from version 3.5 
of Starlab, along with sigma3 and sigma for the automated calculation of cross sections. 



4 TESTS AND COMPARISONS 

To assess the validity of calculations performed with Fewbody, we have compared the results of several scattering experi- 
m ents with the re s ults o f previous studies. For gene ral binary-single interactions, we have compared our results with those 
of iHut fc Bahcaiil l|l983t) : for general binarv-binarv. iMik cglgjjllgggah : and f or detecting close approach distances, we have 



compared with binary-binary calculations performed by 
used extensively and tested thoroughly (see, in particular, 



Bacon et al The scattering facilities in Starlab have been 

Gualandris et alJl2004h . However, there ha s only been one reporte d 



used extensively and tested tnorougniy (see, m particular. Hjualandris et alJlzUU4ll . However, tnere na s only been one reporte c 
comparison between the three-body scattering routine and the A-body routine in the literature dCualandris et al]l2004) 



Below, we perform a new test and show that the two routines agree at a basic level. 



1 See 

2 See 



MtpjT/w^wjmltxdu^^regeau] or search the web for "Fewbody" 
http:// www.manybody.orgl 
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4.1 General binary— single comparison 

iHut fc Bahcaiil lll98 performed one of the most extensive early studies of binary-single star scattering for the equal-mass, 
point-particle case. Figure^Jshows a comparison of the results of 8 x 10^ scattering interactions calculated using Fewbody with 
their Figure 5. Plotted are the total dimensionless cross sections {a/{Tva'^), where a is the binary semimajor axis) for ionization 
(shown by star symbols) and exchange (triangles) as a function of Voo/vc, for the equal-mass, zero eccentricity, point-particle 
case. The dotted lines represent the data from their figure (without error bars), while the straight solid and dashed lines are 
the theoretically predicted cross sect ions for ionizati on and exchange from the same paper. The agreement between the two 
is excellent, although it appears that lHut fc Bahcalll systematically find a slightly larger cross section for ionization. We note, 
however, that the two agree at roughly the one-sigma level. 



4.2 General binary— binary comparison 

The first systematic study of binary-binary scattering was presented bv iMikkolal l)l983ah . He considered binaries w ith equal 
semimajor axes, and stars of equal mass, in the point-particle limit. We have chosen to compare with Table 5 of lMikkolal 
il983ah . which presents sets of scattering experiments performed for several different values of Voo/vc, with impact parameter 
chosen uniformly in area out to the maximum impact parameter found to result in a strong interaction (listed in his Table 3) . 
Only strong interactions were counted, and the eccentricities of the binaries were chosen from a thermal distribution. It should 
be noted, for the sake of completeness, that Mikkola characterised his encounters by their dimensionless energy at infinity. 
Too- The relation between i;oo and Too is Voo/vc = VToo- Mikkola's classification scheme is similar to Fewbody's, the two 
primary difi'erences being: (1) the value of the tidal tolerance, 5, used by Mikkola is 3 x 1 0~^, while the Few body runs use 
S = 10~^; and (2) the criterion use d to assess the dynamical stability of triples is that of iHarringtonl (Il97'jl . a much less 
accurate stability criterion than the bardhn. fc Aarseth i i200lh criterion used by Fewbody. It is therefore expected that the 
classification of Fewbody is more accurate. The binary-binary scattering encounters are classified into five different outcomes. 
The label "undecided" represents an encounter that was deemed to be unfinished after a preset amount of computation time 
- in other words, it could not be classified into one of the four categories of "exchange" , "triple" , "single ionization" , or "full 
ionization". These four outcomes are described in the ricoii = rows of Table |5| Table |21 compares results from Fewbody with 
Mikkola's Table 5. The comparison is also shown graphically in Figure |5] 

Several comments are in order. Looking at the "undecided" column in Table |5] it is clear that Fewbody resolves more 
encounters than Mikkola, yielding roughly half as many undecided encounters. This is a result of both the increased power of 
modern computers - resonant encounters can be integrated longer, and one can use smaller S - and the more accurate triple 
stability criterion available today. In the next column, labelled "exchange" , it is clear that Mikkola finds many more exchange 
encounters than Fewbody. This is thought to be primarily because in this column Mikkola's data include strong interactions, 
which result not only in exchange, but also in preservation. We have not included this type of outcome in the Fewbody results 
because it would have been cumbersome to implement Mikkola's test for a strong interaction. The next column, labelled 
"triples", shows that Mikkola regularly classifies more triples as stable than Fewbody. This results in fewer outcomes labelled 
as "single ionization" , since the test for single ionization occurs after that of triple stability in Mikkola's code. Full ionizations 
can only occur when the total energy of the system is greater than or equal to zero {voo/vc > 1). There is a large discrepancy 
in the number of full ionizations for Voo/vc = 1.225. We are not quite sure of the underlying reason for the discrepancy, but 
think it may be due to the tidal tolerance used, which differs by more than an order of magnitude between the two methods. 
Aside from the systematic discrepancies pointed out above, the two methods agree at a reasonable level, given the differences 
between them. This is especially clear from Figure 13 For all outcomes except full ionization, the methods agree at roughly 
the two-sigma level (the uncertainties shown are one-sigma). 



4.3 Comparison for close-approach distances 

iBacon et alj il996l) presented a more recent and detailed study of binary-binary interactions in the point-particle limit, in 
which close-approach distances were recorded and used to calculate cross sections. In the scattering experiment we have chosen 
for comparison, each binary had equal semimajor axis (oq — ai = a) and zero eccentricity, and all stars had equal mass. 
The impact parameter was chosen uniformly in area out to the maximum impact parameter given b y bmax/a = C/v p o + D , 
where C — 5, and D — 0.6. This expression for the impact parameter is an extension of that used bv lHut fc Bahcal] (Il983ll . 
designed to sample strong interactions adequately. For each encounter, the minimum pairwise close approach distance, rmin, 
was recorded; and from the set, the cumulative cross section calculated. 

FigureElshows a comparison with their Figure 4. T he circles with error bars represent Fewbody data, while the solid-line 
broken power- law is the best fit to the results obtained bv lBacoi ^^lJ ^^^d^ - Ther e is clearly a multiple-sigma discrepancy for 
''min/a^ 0.01. The discrepancy results from the lack of u se bv | B acone ^^Til99(i) of the appropriate algorithm for detecting 
close approach distances with regularization (§ 18.4 of lAarsethf]200^ 7~Sigurdsson has resurrected the original code, and 
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performed a recalculation with smaller timesteps'^. The new result is shown by the dot-dash line. The resulting cross section 
is closer to the Fewbody result, yet still systematically smaller. 

For comparison, we have performed the same calculation using Starlab, shown by the dashed line. The agreement between 
Fewbody and Starlab is excellent. The only discrepancy between the two occurs at rmin/a ~ 1, which represents the weak 
perturbation of binaries due to distant fly-by's. This discrepancy is most likely due to the differing values of the tidal tolerance 
used. For the Fewbody runs, the tidal tolerance was 6 = 10~^, while for the Starlab runs it was 5 = 10^^, causing Starlab to 
numerically integrate some weakly-perturbed binaries that Fewbody treated analytically. The result is a slightly larger Starlab 
cross section for rmin/a ~ 1, as can be seen in the fi gure. 

We should note that the original calculation of iBacon et alJ il996l) was averaged over the range 0.125 < Voo/vc < 0.25, 
while all other results shown in Figure |S] were calculated with Uoo/^c = 0.25. This cannot account for the discrepancy with 
the original calculation, since the inclusion of smaller velocities at infinity will result in more resonant interactions, and hence 
smaller distances of close approach. We have performed calculations with Voo/vc = 0.125 and found that the cross section 
differs from that with Voo/vc = 0.25 by no more than a few per ce nt. 

Finally, we remark that the error in the original calculation of lBacon et al.l |Ii99e1) is only present for small rmin', many of 
the conclusions in their paper are not affected by this error. 



4.4 Comparison between Starlab's three-body and A'^-body scattering routines 

The scattering facilities in Starlab have been used extensively and tested thoroughly ijMcMiUan fc Hutll99fil : fGualandris ei~ai] 
I2OO4I) . However, there is only one reported comparison between sc atters, the three-body scattering routine, and scatter, the 
A'^-body scattering routine, in the literature iCualandris et al]l200'Jl . A simple test, tuned to suit the purposes of this paper, is 
to compare the binaries containing merger products that result from binary-single interactions with those from binary-binary 
interactions designed to mimic binary-single interactions. An obvious choice for the limiting-case binary-binary interaction 
is that in which one binary has an extremely small mass ratio. We performed binary-single runs in which each star had mass 
Mq, radius -Rq, the binary had semimajor axis 1 AU and e — 0, and Voo — lOkm/s. In the binary-binary runs, the binary 
mimicking the single star had a secondary of mass 10~^ ^^Q^ semimajor axis of 20 AU, and e = 0. The results of 10'' runs 
are shown in Figure |7| in which we plot the cumulative fraction of binaries as a function of rp/a, where rp is the pericentre 
distance of the merger binary, and a is the initial binary semimajor axis. The agreement between acatterS (solid line) and 
scatter (dashed line) is good, with both yielding merger binaries with rp strongly concentrated between 0.15 AU and 0.3 AU. 



5 SYSTEMATIC STUDY OF THE COLLISION CROSS SECTION 

To better understand the behavior of the collision cross section, we have systematically studied its dependence on several 
physically relevant parameters. The understanding gained will allow us to reduce the dimensionality of parameter space that 
must be sampled when we later consider MS-star binaries with physically motivated parameters. 



5.1 Dependence on velocity at infinity 

The dimensionless collision cross section {a/{-Ka'^) for binary-single, a/{iT{ao + ai)^) for binary-binary) as a function of 
the relative velocity at infinity, Vao/vc, is shown in Figure |S| for both binary-single interactions (left) and binary-binary 
interactions (right), for several different values of the expansion parameter, /oxp- Circles represent outcomes with one or more 
collisions (two or more stars collide); triangles, two or more (three or more stars collide); and squares, three (four stars collide). 
Red represents runs with /exp = 1; orange, /cxp = 2; green, /oxp = 5; and blue, /cxp = 10. In both experiments (binary-single 
and binary-binary), each star had mass Mq and radius Rq, and each binary had semimajor axis a — 1 AU and eccentricity 
e = 0. The cross section decreases sharply at Vaa/vc = 1, above which resonant scattering is forbidden, and appears to 
approach a constant value, consistent with being purely geometrical. In the resonant scattering regime, below Voo/vc — 1, the 
collision cross section follows the form 1/w^, implying that gravitational focusing is dominant. The ricoii > 1 cross section in 
the resonant scattering regime is quite high, with (j(Doo/«c)^/(7ra^) « 1 for binary-single and cr{voo/vc)^ /{Tv{ao + ii)^) ~ 0.8 
for binary-binary. 

The ncoii > 2 cross section in the binary-single case is about two to three orders of magnitude below that for ricoii > 1, 
depending on /cxp. However, in the binary-binary case, the ricoii > 2 cross section is only down by a factor of a few to 10. 
The reason for the difference is that in the binary-single case, after one collision occurs, there are only two stars left. The two 
remaining stars will either be bound in a binary, or unbound to each other in a hyperbolic orbit. In the case of a bound orbit, 
the two stars are guaranteed to make at least one pericentre passage, and if the merger product in the binary is large enough, 



^ See|http://www.astro.psu.edu/users/steinn/4bod/index.htmI| 
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a collision will occur. In the case of an unbound orbit, the likelihood of a pericentre passage is decreased. In either case, it is 
clear that with only two stars remaining, the complex resonant behavior observed in three- and four-body interactions that 
leads to close approaches will not occur. 

There is a large spread in the ricoii > 3 cross section in binary-binary scattering. This is because it is likely for collision 
products to suffer subsequent collisions given their increased size, implying that the ricoii > 3 cross section should vary as 
/cxp- The ricoii > 3 cross section varies from a factor of a few to two orders of magnitude below that for ricoii > 2. 

Finally, we note that the spread in the Wcoii > 2 binary-binary cross section is a factor of about four, essentially 
independent of v^o for Voo/vc < 1, as /cxp varies over an order of magnitude. The cross section is therefore not a particularly 
sensitive function of the unknown expansion parameter /exp, and, if it is valid to parametrize the size of collision products in 
this simplified manner, implies that our results for the properties of merger populations are relatively robust. 

5.2 Dependence on the ratio of stellar radius to binary semimajor axis 

The collision cross section varies as for Voo/vc < 1, the regime relevant to interactions involving hard binaries in the 

cores of globular clusters. Therefore, we can choose a single value for Voo when exploring the dependence of the collision 
cross section on other physically relevant parameters, thereby reducing the dimensionality of parameter space that must be 
sampled. For the remainder of this section, we set Voo/vc = 0.1, which corresponds to typical binary-single and binary-binary 
interactions involving hard binaries in a globular cluster core, with Voo ~ lOkm/s, stars of mass Mq, radius Rq, and binaries 
with a = 0.1 AU. 

Figure |n| shows the normalized, dimensionless collision cross section, a{voc / / {-ko?) for binary-single scattering (left), 
<j{voo/vcY /{^{ao + O'lf') for binary-binary scattering (right), as a function of the ratio of stellar radius to binary semimajor 
axis, R/a, for different values of the expansion parameter, /oxp- Circles represent outcomes with one or more collisions; 
triangles, two or more; and squares, three or more. Red represents runs with /exp = 1; orange, /cxp = 2; green, /cxp = 5; and 
blue, /exp = 10. In both experiments (binary-single and binary-binary), each star had mass Mq and radius R, each binary 
had semimajor axis a = 1 AU and eccentricity e = 0, and the relative velocity at infinity was set to Woo/fc = 0.1. Calculations 
were performed down to R/a = 10~^ - which corresponds to the extreme case of binaries with semimajor axis 10 AU composed 
of black holes of mass Mq - but no collisions were found below R/a « 10~®. For ricoii > 1, the calculation corresponds to the 
simpler task of recording minimum close approach distances, as can be seen by comparing the binary-binary panel (right) to 
Figure |S| The ricoii > 2 and ncoii > 3 collision cross sections decrease more sharply than the ricoii > 1 cross section as R/a 
decreases. 

It is clear that multiple collisions are unlikely for R/a< 0.001, which corresponds roughly to stars of radius Rq in binaries 
with semimajor axis 1 AU. We therefore expect that multiple collisions in binary interactions are relevant only for MS stars 
in binaries tighter than ~ lAU, white dwarfs in binaries tighter than ~ IRq, and neutron stars in binaries tighter than 
~ 10"' km. We caution that relativistic effects may need to be included when considering close approaches of neutron stars. 
However, the limits quoted should serve as a rough guide. 

We have held the stellar masses fixed at Mq, while varying their radii over a large range. For MS stars, it is more realistic 
to adopt a reasonable mass-radius relationship, which we do in § El for several sets of masses. 

5.3 Dependence on mass ratio 

In binary intera ctions involving star s of different masses, there is a strong tendency for the lightest star(s) to be ejected 
quickly (see, e.g.. lHeggie fc Hxill2003h . One would expect, then, that resonant behavior, and the likelihood of collisions, would 
be decreased when one or more of the stars involved are light. To test this prediction, we have calculated the collision cross 
section during binary-single and binary-binary scattering for a range of mass ratios. In both experiments, each binary had 
one star with mass Mq and the other with mass qMq. For the binary-single case, the incoming single star had mass Mq. 
Each star had radius Rq, each binary had semimajor axis a = 1 AU and eccentricity e = 0, and the relative velocity at infinity 
was set to Voc/vc = 0.1. We normalize the cross section, as usual, by multiplying by (vao/vcf^ , and, in doing so, inadvertently 
introduce a dependence on the mass ratio, q, in v^. To remove it, we also multiply by a function of q alone that has the same 
dependence on q as v1, from eqs. and Q, and is normalized to 1 at g = 1. For binary-single interactions this function is 
2q{2 -\- g)/(3(l -f q)); for binary-binary interactions it is 2g/(l + q). The collision cross sections are shown in Figure [101 for 
binary-single (left) and binary-binary (right), as a function of q, for different values of the expansion parameter, /cxp- Circles 
represent outcomes with one or more collisions; triangles, two or more; and squares, three or more. Red represents runs with 
/cxp = 1; orange, /exp = 2; green, /cxp = 5; and blue, /cxp = 10. As expected, the collision cross section is smaller for q < 1. 
However, it decreases quickly, and for g<0.1 becomes approximately constant, implying that the test particle limit has been 
reached. What is most striking is that the collision cross section is decreased by no more than a factor of a few for small q, 
despite the tendency for lighter stars to be ejected quickly. It should be noted that in this experiment we have kept the radii 
of all stars fixed at Rq. 
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6 RESULTS FOR TYPICAL BINARIES 

We now turn from a slicing of parameter space to a discrete sampling, by considering binaries with sets of parameters typical 
of those found in the cores of globular clusters. We first present results for binary-single interactions, and then binary-binary. 

6.1 Binary— single scattering experiments 

We consider only MS stars with masses 0.5 Mq, 1.0 Mq, or 1.2 Mq. We adopt the mass-radius relationship R = Rq{M/Mq), 
which is a reasonable approximation for MS stars of mass ~ 1 Mq . We study five different mass combinations, labelled A 
through E, with a range of semimajor axes, 0.05 AU < a < 3.0 AU, for each. In all cases we use = lOkm/s. This choice of 
parameters covers a range of binary binding energies from ~ 1 kT (the hard-soft boundary) in a typical globular cluster core, 
to ~ 10^ kT, corresponding to a close binary (a ~ IO-Rq). The thermal energy kT is defined by the relation ^kT — ^{m)a^, 
where (m) is the average stellar mass, and a is the one-dimensional velocity dispersion. The details of each run are presented 
in Table |1J including run name; the number of scattering interactions performed. A''; the masses of the binary members, moo 
and moi; the mass of the intruder, mi; the binary semimajor axis, a; and the ncoii > 1 cross section. 

In order to study the dependence of the collision cross section on the expansion parameter, /cxp, without performing 
calculations for each value of /cxp considered, we have adopted an approach that allows us to calculate multiple collision cross 
sections for any value of /oxp based on the results of calculations for one value of /cxp- We set /cxp ~ 1, and consider the 
properties of merger binaries formed. A binary containing a merger product will be a triple-star merger if the pericentre of 
the binary, rp, is approximately less than the radius of the collision product. Rep = fc^p{Ri + R2), where Ri and R2 are the 
radii of the two stars that merged to form the collision product. First we calculate A'coii, the total number of outcomes that 
resulted in either merger binaries or triple mergers with /cxp ~ 1. We then calculate A^scoii, the number of triple mergers, for 
a different value of /cxp, as the number of triple mergers for /cxp ~ 1, plus the number of merger binaries with Vp < Rep. 
Defining fr = iVscoii /A'coii, the triple-star merger (wcoii > 2) cross section for /cxp is simply arifc^p) = /Tcrcoii(/exp = 1). 

Some remarks about this approach are in order. We ignore merger escapes, and argue that an outcome labelled as a 
merger escape is unlikely to become a triple merger even if the first merger product expands. Before it escapes, the third 
star can approach the expanded merger at most once, and, if it does, it is likely to have a sufficiently high speed at close 
approach to fully traverse the tenuous envelope of the expanded merger product. On the contrary, in a merger binary, even if 
the third star initially has a high pericentric speed, it will eventually be captured through gradual energy loss after repeated 
traversal. Of course, an escaping third star may lose sufficient energy after traversal so that the entire system becomes bound, 
and eventually be captured. A more precise treatment would be to run calculations for each value of /exp, but, as mentioned 
above, we are adopting the simpler, less computationally expensive approach here. When two MS stars collide and their merger 
product expands, the resulting object does not possess a well-defined boundary and, in general, is not spherically symmetric; 
/cxp is thus an effective, averaged quantity, which serves well enough the purpose of our first study. 

6.1.1 Collision cross sections 

The ricoii > 1 cross sections are listed in the last column of Table |1| The cross sections from runs A, B, and C are also 
shown as a function of the initial binary semimajor axis, a, in Figure fTTI In the range of MS masses of interest for globular 
clusters, the collision cross sections show only a weak dependence on masses, slightly more pronounced at small a. The cross 
section increases from case A to C as the mass ratios of the stars decrease, due to the dependence of the normalized cross 
section on Voc/vc and hence on the mass ratio. For a hard binary with a ~ 1 AU, the normalized collision cross section 
is comparable to the geometric cross section of the initial binary (i.e., a-con{vac/vc)'^ ~ iva^). This is because most strong 
interactions are resonant, and most resonances lead to at least one collision. For a<0.1AU, the collision cross section can 
be up to an order of magnitude greater than the geometric cross section. Indeed, for very small values of a, even a small 
perturbation of a highly eccentric orbit by a distant encounter can induce a binary merger. About 20 to 35% of the initial 
binaries with a — 0.05 AU in Table |1| have pericentre distances less than 3Rq. Our results for these very tight binaries 
are therefore somewhat artificial, since in reality tidal circularization effects are likely to modify the distribution of initial 
eccentricities, and our simple assumption of a thermal initial distribution is no longer justified. 

6.1.2 Properties of the merger binaries 

Of particular interest are binary-single interactions that result in binaries containing merger products. The distributions of 
their properties are relevant to observations of BSs in the cores of globular clusters. Figures 1121 and 1131 show the orbital 
parameters of the merger binaries produced in the two representative runs A300, for a wide initial binary, and B005, for 
a very tight initial binary. The envelope of the distribution follows curves of constant angular momentum, consistent with 
angular momentum conservation during the interaction. The total angular momentum of the system is the sum of the initial 



© 2004 RAS, MNRAS 000, 000-000 



12 J. M. Fregeau, P. Cheung, S. F. Portegies Zwart, & F. A. Rasio 



internal angular momentum of the binary and the initial angular momentum of the binary-single hyperbolic orbit, added 
vectorially. The spread in angular momentum spanned by the distributions is due to averaging over the relative orientation 
of the two separate angular momenta, the range of initial eccentricities of the binary, and range of impact parameters 
used. Curves of constant angular momentum are plotted in Figure fT^ for the values J/ Jo = 0.2, 0.5, 1.0, and 2.0, where 
Jo = nbvoo + i-Lb[GA4ba{l — e^)]^^'^ is the angular momentum of the system such that the pericentre distance of the initial 
hyperbolic orbit is 1.0 AU (i.e., b = rp(l + 2GM /rpV^)^^^ with rp = 1.0 AU). (Here fj, and M are the reduced and total 
mass of the binary-single system, and fib and M(, are the reduced and total mass of the binary.) The vertical dashed line in 
Figures fT^ and [T^ is the hard-soft boundary with respect to field stars of mass 1 AJq with one- dimensional velocity dispersion 
lOkm/s. Histograms of final semi-major axes and eccentricities are shown in Figures mi and [T51 The dotted lines in Figure 
1151 are properly normalized thermal eccentricity distributions. 

Typically, more than 90% of the merger binaries have final semimajor axis, a', larger than initial, a. On average, a' /a « 5. 
While most remain hard binaries, a small fraction become soft, with a few having a' as large as ~ 100— 1000 AU. This softening 
comes from the somewhat counter-intuitive result that collisions produce, on average, an increase in the orbital energy of the 
system (while the total energy, including the binding energy of the collision product, is of course conserved on a dynamical 
time-scale, i.e., until some of the internal energy released through shocks can be radiated away by the fluid). To illustrate 
this, consider a trivial example in which two identical stars of mass m are released from rest at some distance r and collide 
head-on, forming a stationary merger product at the centre of mass. The orbital energy of the system increased by Gm^/r in 
the process. More relevant to our results, but still somewhat artificial, consider an initial binary with a very high eccentricity, 
so that the two members almost collide at pericentre. A small perturbation through a distant encounter can induce a merger 
of the binary (implying that its orbital binding energy disappears), while only weakly affecting the orbit of the perturber. 

The eccentricity distributions of merger binaries always remain close to thermal, although a slight excess of highly eccentric 
orbits is seen for wider initial separations (compare run A300, with a — 3AU, and B005, with a = 0.05 AU, in Figure [T^ . 
The average value of e' for runs A300 and B005 is 0.77 and 0.68, respectively, while that of a thermal distribution is 2/3. 
It is interesting to note that other calcul ations of small- systems have yie lded binaries with an excess of high eccentricity 
systems in a nearly thermal distribution ijPortegies Zwart fc McMillanll200(l l. 



6.1.3 Three-star mergers 

Three-star mergers happen primarily when the pericentre distance of a merger binary is approximately smaller than the radius 
of the merger remnant. Cumulative distributions of pericentre distances from all A runs are shown in Figure \Wi For radii of 
first collision products in the range ~ 5-10 -Rq (/cxp ~ 2.5-5), we find triple collision fractions anywhere from a few percent 
up to 50%, depending strongly on the initial binary semi-major axis a. Clearly, triple collisions occur often, particularly during 
encounters with very hard binaries. If we consider the later expansion of the collision product on the giant branch (with radius 
up to > 100 AU), a triple collision becomes almost inevitable, except for only the widest initial binaries. 

Denoting the value of -Rep at which /t = / by f?/, we determine the critical radii of merger products corresponding to 
a given triple collision fraction, Ro.os, Ro.i, R0.5, -R0.9 and iio.95, using simple linear interpolation. These are plotted as a 
function of the initial binary separation a in Figure IT7I The error bars in Figure IT7I are estimated by dividing the uncertainty 
in /t by the slope of the /t vs Rep curve at Rep — Rf, i.e.. 



- {df/dRep)n, ■ 

We see that all the lines in Figure 1171 are nearly parallel and with a slope close to unity. The same holds true for mass 
combinations B through E as well. Thus we have approximately Rf oc a and the relationship can be specified by a single 
quantity Rf/a for each value of /t. These have been estimated using a least-squares fit with weights inversely proportional 
to the size of the error bars. Since hydrodynamic calculations have shown that _Rcp is unlikely to be larger than ~ 30 times 
the original stellar radius, according to Figure IT7I the most relevant range corresponds to /t^O.5 (although the full range 
up to /t ~ 1 will be relevant if the later expansion of the merger product on the red giant branch is considered). In Figure 
ITslwe plot Rf/ a as a function of /t. It is clear that Rf/a is directly proportional to /t over the range of interest. For run 
A (equal mass case), the proportionality constant is 1.61 ± 0.01. Consequently, the relation between -Rep, fr and a for this 
particular mass combination may be written 

-Rep ~ 1.6a/T , (9) 
where 0.05 AU < a < 3.0 AU. Turning to different mass combinations we find results similar to eq. @, and so can write 
-Rep = CafT , (10) 
where C depends only on the stellar masses. Table |S| shows C for the five mass combinations we have explored. 
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6.2 Binary— binary scattering experiments 

For the sake of convenience, we use the abbreviations listed in Table |5| to refer to certain binary-binary outcomes. In the 
abbreviated form, the letters S, D, T, and Q denote a single star, double-star merger, triple-star merger, and quadruple- star 
merger, respectively, and we have also chosen to use parentheses instead of square brackets. Each run we do involves MS stars 
of either 0.5 Mq or 1.0 Mq and binary semimajor axes of either 1.0 AU or 0.1 AU. In all runs we set Voa ~ 10 km/s, as in the 
binary-single case. The properties of each run are listed in Table including the mass of each star, rrnj, the semimajor axis 
of each binary, a, and the normalized cross sections for strong interactions and at least one collision to occur. 

To study the dependence of the outcomes on the expansion parameter, /cxp, we have performed separate calculations 
for each value of /oxp considered. For binary-binary interactions, the dynamics do not reduce to the trivial analytical case of 
two-body motion after one collision has occurred, and so it is not possible to use the simple approach of tracking pericentre 
distances in merger binaries as we did for the binary-single case. It should be noted, however, that we apply the simple 
expansion factor prescription for the radius of a merger product, -Rmorgor = .fcxp{Ri + R2), where Ri and R2 are the radii 
of the merging stars, to every merger, regardless of whether the merging stars are unperturbed MS stars or merger products 
themselves. The simplicity of this prescription allows us to study the dependence of our results on only one parameter, /cxp, 
which can thus be considered an effective expansion parameter, averaged over all types of mergers. A more realistic approach 
that adopts separate expansion parameters for different types of mergers is feasible, but beyond the scope of this study. 

6.2.1 Collision cross sections 

The normalized cross sections for strong interactions, (Jstrong, and for at least one collision, ctcoii, for our binary-binary runs 
are listed in the last two columns of Table |S| A strong interaction is defined to be one in which the final configuration is 
different from the initial configuration (i.e., anything but preservat ion) , or a preservation resulting from a resonant encounter. 
The test for a resonant encounter is that of lHut fc Bahcalll l)l98.'^ ). wherein the mean square distance between pairs of stars 
is checked for multiple minima. 

Comparing the results from run I (ao = ai = 1 AU) with run II (ao = ai =0.1 AU), we see that (Jcoii is a larger fraction of 
""strong for run II, consistent with our findings in i; l5.2l for small R/a. Comparing run II with run IV, we see that introducing a 
non- unity mass ratio does not seem to affect (Jstrong, but slightly lowers ctcoii. By calculating the branching ratio for outcome 
X involving collisions - defined as fx — Nx/NcoU where A'coii is the number of outcomes that result in collisions and Nx 
is the number of those that result in outcome X - the value of CTcoII for a particular run can be used to calculate the cross 
section for outcome X, according to the simple relation ax = /xctcoII- 

6.2.2 Properties of merger products 

In Figures^] and ism we show the branching ratios for several outcomes as a function of /oxp. That is, we plot the fraction of 
outcomes involving at least one collision that result in various configurations containing double-star, triple-star, and quadruple- 
star mergers. Figure [T^ shows results from run I (ao = ai = 1 AU) and Figure[5U|shows results from run II (ao — ai —0.1 AU). 
The upper left panel in each shows the branching ratios for outcomes of two unbound double-star mergers, labelled DD, and 
two double-star mergers in a binary, labelled (DD); the upper right, a quadruple-star merger, labelled Q; the lower right, a 
triple-star merger bound to the remaining single star, labelled (TS); and the lower left, the combined branching ratio for any 
outcome involving a merger of three or more stars, labelled T/Q. 

From Figure IT^ we see that, even for encounters involving wider binaries, the branching ratio for more than two stars to 
merge is significant - as high as ~ 5%. When one considers tighter binaries, as in Figure 1^ the branching ratio increases to 
~ 40%. The dependence on initial semimajor axis is as expected - all branching ratios for mergers are increased in run II over 
run I. The dependence on /oxp is also as expected. As the expansion factor is increased, more multiple mergers occur, leading 
to an increase in the branching ratios for triple-star and quadruple-star mergers, and a decrease in those for double-star 
mergers. 

The distributions of orbital parameters for all four types of binaries, (DS)S, D(SS), (DD), and (TS) are plotted in Figures 
EH and El for runs I and II, with /exp = 5. From these figures, we see that (DS) binaries form with semimajor axes comparable 
to, and only slightly greater than, the semimajor axes of their progenitor binaries (except for the case (DS)S, where a' can 
be significantly larger than a for large e'), and with an eccentricity distribution that does not appear to be inconsistent with 
thermal. The data are more sparse for the (DD) and (TS) cases, but their orbital parameters appear to be comparable to 
those of (DS) binaries. 

The three outcomes TS, (TS), and Q (labelled T/Q collectively in Figures and ESJ, are responsible for the production 
of BSs of mass > 2Mq in our runs. The branching ratio for T/Q appears to increase almost linearly with /cxp in the range 
considered, for all runs performed. Linear fits for the branching ratio of T/Q as a function of /cxp (obtained by least squares 
fitting) are provided in Table |7| 
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7 SUMMARY, CONCLUSIONS, AND FUTURE DIRECTIONS 

We have performed several sets of binary-single and binary-binary scattering experiments, and studied the likelihood of 
(multiple) collisions. We have presented collision cross sections, branching ratios, and sample distributions of the parameters 
of outcome products. Results reported in this paper, particularly cross sections, may be employed in both analytical and 
numerical calculations. 

In the gravitational focusing regime, relevant to hard binaries in globular cluster cores, the likelihood of collisions during 
binary interactions is quite high. For solar mass main-sequence (MS) stars in 1 AU binaries, the normalized cross section 
for at least one collision to occur during a binary-single or binary-binary interaction (ricoii > 1) is essentially unity, with 
o'{vcx,/vc)'^ /{TTa^) ~ 1 for binary-single and o"(woo/'Jc)^/(7i"(ao + ai)^) ~ 1 for binary-binary. The collision cross section depends 
strongly on the ratio of stellar radius to binary semimajor axis, but is reasonably high even for MS stars of approximately 
solar mass in orbits of ~ 1 AU. Perhaps counter to intuition, the collision cross section is not particularly sensitive to binary 
mass ratio, dropping by only a factor of a few in the test-particle limit when the stellar radii are kept fixed. We also found 
that the multiple collision (ricoii > 2) cross section is quite high, only a factor of ~ 10 lower than the ncoii > 1 cross section for 
binary-binary interactions. It is also not a particularly sensitive function of the expansion parameter, /exp, varying by a factor 
of a few as /cxp is varied by an order of magnitude. This implies that studies using this one-parameter model for the radius 
of a collision product are reasonably robust in spite of the large uncertainties in the physics. For typical binaries in globular 
cluster cores, we have shown that collisions of more than two stars during binary-single and binary-binary interactions are 
likely, with branching ratios for triple-star mergers of ~ 5% for binary-single and ~ 10% for binary-binary. 

We have introduced Fewbody, a new numerical toolkit for simulating small- A'' gravitational dynamics that is particularly 
suited to performing scattering interactions. We have shown that it produces results in good agreement with several previous 
numerical studies of binary-single and binary-binary scattering, as well as with the Starlab software suite. Instead of using 
cross sections and simple recipes for binary interactions in globular cluster evolution codes, one may use Fewbody to perfo rm 
them directly. We have adopted this approach with our Monte Carlo globular cluster evolution code ijFregeau et ahlEoO^ll . 

It is clear from our results that collisions of more than two stars during binary interactions are a viable pathway for 
creating blue straggl ers (BSs) with masses more than twice the MS turnoff mass, such as those observed in NGC 6397 
JSepinskv et al.l2000l) . These massive BSs may also be formed via recycling - in other words, a binary containing a BS may be 
formed via a binary interaction, and the BS may later merge with another star in a subsequent binary interaction, creating a 
more massive BS. We are in the process of creating a more detailed model, based on a Monte Carlo binary population study, 
which incorporates b oth channels to study the formation of massive BSs. Such studies are needed to help interpret current 
BS observations (see ISills fc Bailvnlll999l : ISills et aklboool) . and the large databases of BS properties, including many new 
spectroscopic mass measurements, that will soon be available (M. Shara, private communication). 

The expansion of merger products has been treated here in a simplified manner, using a single expansion parameter, /cxp. 
As observations of BSs become more detailed and more numerous, including details of their internal properties, it becomes 
necessary to treat collisions in a more accurate way. Full smoothed particle hydrodynamics (SPH) calculations are quite 
computationally prohibitive, taking up to several hours to perform a single merger. However, there are faster, approximate 
approaches that capture the essential physics of the hydrodynamic merger process. One such approach is the fluid-sorting 
algorithm, which utilises th e property that the fiuid in merger products must rearrange itself according to specific entropy 
iLombardi et al.lll995l . E'o02^ . The Make Me A Sta r (MMAS) software developed by L ombardi and collaborators implements 
this procedure, and is freely available on the web l|Lombardi et alJIlQQfil . 120021 . 1200.^ 1 . We have begun to replace the simple 
merger module in Fewbody with a call to MMAS. The result should be much more accurate predictions for the properties of 
(multiple) merger products. 
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Table 1. Possible outcomes of binary-single star encounters, ordered by the number of collisions, n^all- Brackets enclose two objects 
which are bound to each other, while colons represent physical collisions. For simplicity, we have only listed the outcomes that would 
result from indistinguishable stars. 



"coll 


symbol 


description 





[. .] . 


preservation or exchange 





• • • 


ionization 





[[• •] •] 


stable hierarchical triple 


1 


[.:. .] 


binary containing a two-star merger 


1 


• :• • 


two-star merger and single star 


2 


• 


three-star merger 



Table 2. Possible outcomes of binary-binary star encounters, ordered by the number of collisions, n^^ll- Brackets enclose two objects 
which are bound to each other, while colons represent physical collisions. For simplicity, we have only listed the outcomes that would 
result from indistinguishable stars. Listed in the third column are the abbreviations used in the paper to refer to various outcomes. 



"coll 


symbol 


abbreviation 


description 





[• •] [• •] 




preservation or exchange 









single ionization 





• • • • 




full ionization 





[[• •] •] • 




stable hierarchical triple and single star 





[[[• •] •] •] 




stable hierarchical quadruple 





[[• •] [• •]] 




stable quadruple composed of two binaries 


1 


[• •] •:• 


(SS)D 


binary and two-star merger 


1 


[•:• •] • 


(DS)S 


single star and binary containing two-star merger 


1 


• :• • • 


DSS 


two-star merger and two single stars 


1 


[[.:. .] .] 


({DS)S) 


stable hierarchical triple with two-star merger in inner binary 


1 


[[. .] .:.] 


({SS)D) 


stable hierarchical triple with two-star merger in outer binary 
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[•:• •:•] 


(DD) 


binary composed of two two-star mergers 


2 


[•:•:• •] 


(TS) 


binary containing a three-star merger 


2 


• :• •:• 


DD 


two two-star mergers 


2 


• • 


TS 


three-star merger and single star 


3 




Q 


four-star merger 
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Figure 1. Schematic representation of the binary-tree algorithm used by Fewbody. A circle containing a number represents a star. The 
set of binary trees is shown flattened here, as it is before processing, so that each star is the top-level node of a one-node tree. 




Figure 2. Schematic representation of the binary-tree algorithm used by Fewbody. A circle containing a number represents a star, while 
an empty circle represents a general parent node. The set of binary trees is shown after the first stage of processing, with stars 1 and 2 
replaced by their parent node, which contains the dynamical information pertaining to the centre of mass of the 1-2 binary, as well as 
all phase and orientation information. For classification, the replacement of stars 1 and 2 by their parent node simply means that they 
are bound to each other with the smallest semimajor axis. For hierarchy isolation, it would also mean that the 1-2 binary is sufficiently 
weakly perturbed that it can be treated analytically, and is stable in the sense that its two members will not collide at pericentre. 
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Figure 3. Schematic representation of the binary-tree algorithm used by Fewbody. A circle containing a number represents a star, while 
an empty circle represents a general parent node. The set of binary trees is shown after the second stage of processing, with the 1-2 
centre of mass and star 3 replaced by their parent node. For hierarchy isolation, this replacement is quite rare, as it would require that 
the triple be not only dynamically stable, but also sufRciently hierarchical that its evolution could be treated analytically. 




Voo/Vc 



Figure 4. Comparison of Fewbody with Figure 5 of lHut &: BahcalU ^^sj): total cross sections for binary-single scattering for the equal- 
mass, zero eccentricity, point-particle case. A total of 8 X 10^ scattering experiments were used to create this figure. The dotted lines 
represent the data from lHut fc Bahcalll il983l) . while the straight solid and dashed lines are the theoretically predicted cross sections for 
ionization and exchange from the same paper. Data points are from Fewbody. The agreement between the two is excellent. 

Table 3. Comparison of Fewbody with Table 5 of lMikkolal il983aft : fraction of strong binary-binary interactions that result in various 
outcomes. In each binary— binary interaction the stars had equal masses and were assumed to be point particles, the binaries had 
equal semimajor axes, and the eccentricities were drawn from a thermal distribution. The data are normalized to 100 total scattering 
experiments. The results are also shown graphically in Figure IBl 
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method 
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exchange 


triple 


single ionization 


full ionization 


total 


0.316 


Mikkola 


4.7 ±1.2 


11.0 ± 1.9 


24.3 ±2.8 


60.0 ±4.5 


0.0 ±0.0 


300 




Fewbody 


1.6 ±0.2 


6.0 ±0.3 


22.8 ±0.7 


69.7 ± 1.2 


0.0 ±0.0 


5225 


0.500 


Mikkola 


2.7 ±0.9 


7.3 ± 1.6 


20.3 ±2.6 


69.7 ±4.8 


0.0 ±0.0 


300 




Fewbody 


1.4 ±0.2 


6.7 ±0.4 


17.2 ±0.6 


74.7 ±1.3 


0.0 ±0.0 
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0.707 


Mikkola 


0.7 ±0.5 


9.3± 1.8 


11.7 ±2.0 


78.3 ±5.1 


0.0 ±0.0 


300 




Fewbody 


0.5 ±0.1 


7.7 ±0.5 


6.8 ±0.5 


85.0 ± 1.6 


0.0 ±0.0 


3303 


0.866 


Mikkola 


0.7 ±0.5 


16.7 ± 2.4 


5.0 ±1.3 


77.7 ±5.1 


0.0 ±0.0 


300 




Fewbody 


0.1 ±0.1 


8.2 ±0.5 


2.3 ±0.2 


89.4 ± 1.5 


0.0 ±0.0 


3827 


1.000 


Mikkola 


0.0 ±0.0 


9.3 ± 1.8 


1.7 ±0.7 


89.0 ±5.4 


0.0 ±0.0 


300 




Fewbody 


0.1 ±0.0 


6.4 ±0.4 


0.6 ±0.1 


92.7 ± 1.6 


0.2 ±0.1 


3499 


1.225 


Mikkola 


0.0 ±0.0 


8.3± 1.7 


0.7 ±0.5 


73.7 ±5.0 


17.3 ±2.4 


300 




Fewbody 


0.0 ±0.0 


4.6 ±0.3 


0.2 ±0.1 


92.1 ± 1.5 


3.1 ±0.3 


3969 
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Figure 5. Comparison of Fewbody (solid lines) with Table 5 of lMikkoIal ll983a^ (dotted lines): fraction of strong binary— binary interactions 
that result in various outcomes. In each binary-binary interaction the stars had equal masses and were assumed to be point particles, 
the binaries had equal semimajor axes, and the eccentricities were drawn from a thermal distribution. Circles represent outcomes that 
were undecided after a preset maximum computation time, squares represent exchanges, diamonds represent stable hierarchical triples, 
upward-pointing triangles represent outcomes that resulted in one binary being disrupted, and downward-pointing triangles represent 
outc omes that resulte d in both binaries being disrupted. The solid lines represent Fewbody data, while the dotted lines represent data 
from lMikkoial EUs^)- The results are also presented in Table 1^ 




^min/^ 

Figure 6. Comparison of Fewbody with Figure 4 of lBacon et alJ ^^^S)- cumulative cross section for the distance of closest approach 
in binary— binary scattering for the equal-mass, zero-eccentricity, equal-semimajor-axis case. The stars are assumed to be point particles 
an d Voo/vc = 0.25. A total of 1.5 X 10* scattering experiments were used to create this figure. The broken power-law is the best fit given 
bv lBacon et al.l il99ff) to their original results, while the dot-dash curve is Si gurdsson's recalculat ion. The dashed curve shows the results 
obtained using Starlab. There is a clear discrepancy between Fewbody and lBacon et all Jl996l) . and even the recalculation. However, 
Fewbody and Starlab agree quite well. 
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Figure 7. Comparison between scatters, Starlab's three-body scattering routine (solid line), and scatter, its A^-body scattering routine 
(dashed line). Plotted is the cumulative fraction of binaries as a function of rp/a, where rp is the pericentre distance of the merger binary, 
and a is the initial binary semimajor axis. For the binary-single runs, each star had mass Mq, radius Rq, the binary had semimajor 
axis lAU and e = 0, and Voo = lOkm/s. In the binary-binary runs, the binary representing the single star had a secondary of mass 
10~^ semimajor axis of 20 AU, and e = 0. The agreement between the two methods is excellent, and in either case, rp is strongly 

concentrated between 0.15AU and 0.3 AU. 
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Figure 8. Cross section for physical collisions in binary-single (left) and binary-binary (right) scattering as a function of the relative 
velocity at infinity, for different values of the expansion parameter, /exp - Circles represent outcomes with one or more collisions; triangles. 



two or more; and squares, three or more. Red represents runs with /exp = 1; orange, /e; 



2; green, /exp = 5; and blue, /exp = 10. 
In both experiments (binary— single and binary— binary ) , each star had mass Mq and radius Rq, and each binary had semimajor axis 
a = 1 AU and eccentricity e = 0. The cross section decreases sharply at the critical velocity, tic, above which resonant scattering is 
forbidden. 
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Figure 9. Normalized cross section for physical collisions in binary-single (loft) and binary— binary (right) scattering as a function of the 
ratio of each star's radius to each binary's semimajor axis, R/a, for different values of the expansion parameter, /cxp- Circles represent 
outcomes with one or more collisions; triangles, two or more; and squares, three or more. Red represents runs with /oxp = 1; orange, 



2; green, /e; 



5; and blue, /exp 



radius R, each binary had semimajor axis a 
Calculations were performed down to R/a = 



10. In both experiments (binary— single and binary-binary), each star had mass Mq and 
= 1 AU and eccentricity e = 0, and the relative velocity at infinity was set to Vao/^c = 0.1. 
10~^, but no collisions were found below R/a ^ 10~^. 
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Figure 10. Normalized cross section for physical collisions in binary-single (left) and binary-binary (right) scattering as a function of 
mass ratio, q, for different values of the expansion parameter, /exp. Circles represent outcomes with one or more collisions; triangles. 



two or more; and squares, three or more. Red represents runs with fe- 



1; orange, 



: 2; green, /o; 



5; and blue, fc: 



10. In 



both experiments (binary-single and binary-binary), each binary had one star with mass Mq and the other with mass qMq. For the 
binary— single case, the incoming single star had mass Mq. Each star had radius Rq, each binary had semimajor axis a = lAU and 
eccentricity e = 0, and the relative velocity at infinity was set to v^c/vc = 0.1. 
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Table 4. Parameters of the binary-single runs, including the number of scattering interactions performed, A'^; the masses of the binary 
members, moo and moi; the mass of the intruder, mi; the binary semimajor axis, a; and the ricoU > 1 cross section. 



run 


N 


"^00 (-'^^o) 


moi (Mq) 


mi (A/q) 


a(AU) 


/ \ 2 

■KO^ \ Vc / 


A005 


15054 


1.0 


1.0 


1.0 


0.05 


6.4 ± 0.1 


AOlO 


30228 


1.0 


1.0 


1.0 


0.1 


5.72 ± 0.07 


A020 


15222 


1.0 


1.0 


1.0 


0.2 


4.41 ± 0.08 


A050 


18158 


1.0 


1.0 


1.0 


0.5 


2.87 ± 0.07 


AlOO 


37625 


1.0 


1.0 


1.0 


1.0 


1.94 ±0.04 


A300 


21427 


1.0 


1.0 


1.0 


3.0 


0.68 ± 0.04 


B005 


17619 


1.0 


0.5 


1.0 


0.05 


9.3 ±0.2 


BOlO 


30408 


1.0 


0.5 


1.0 


0.1 


7.3 ± 0.1 


B020 


17969 


1.0 


0.5 


1.0 


0.2 


5.4 ±0.1 


B050 


18676 


1.0 


0.5 


1.0 


0.5 


3.0 ±0.1 


BlOO 


39739 


1.0 


0.5 


1.0 


1.0 


1.62 ± 0.05 


B300 


28544 


1.0 


0.5 


1.0 


3.0 


0.54 ± 0.05 


COOS 


17696 


0.5 


0.5 


1.0 


0.05 


12 ±2 


COlO 


35791 


0.5 


0.5 


1.0 


0.1 


9.3 ± 0.1 


C020 


18284 


0.5 


0.5 


1.0 


0.2 


6.6 ± 0.2 


C050 


19467 


0.5 


0.5 


1.0 


0.5 


3.3 ± 0.1 


ClOO 


49032 


0.5 


0.5 


1.0 


1.0 


1.96 ± 0.07 


C300 


33464 


0.5 


0.5 


1.0 


3.0 


0.58 ±0.07 


D005 


12530 


1.0 


1.0 


0.5 


0.05 


2.5 ±0.1 


DOlO 


12555 


1.0 


1.0 


0.5 


0.1 


2.2 ±0.1 


D020 


12610 


1.0 


1.0 


0.5 


0.2 


1.75 ±0.1 


D050 


12780 


1.0 


1.0 


0.5 


0.5 


1.16 ± 0.08 


DlOO 


15672 


1.0 


1.0 


0.5 


1.0 


0.73 ± 0.07 


D300 


14185 


1.0 


1.0 


0.5 


3.0 


0.26 ± 0.04 


E005 


60252 


1.0 


1.0 


1.2 


0.05 


8.0 ±0.2 


EOlO 


60504 


1.0 


1.0 


1.2 


0.1 


6.9 ±0.2 


E020 


61008 


1.0 


1.0 


1.2 


0.2 


5.4 ±0.2 


E050 


72947 


1.0 


1.0 


1.2 


0.5 


3.5 ±0.1 


ElOO 


75894 


1.0 


1.0 


1.2 


1.0 


2.2 ± 0.1 


E300 


100200 


1.0 


1.0 


1.2 


3.0 


0.80 ± 0.05 




a (AU) 

Figure 11. The normalized n^^n > 1 cross section as a function of initial semimajor axis for runs A, B and C. 
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Figure 12. Distribution of the scmimajor axes and eceentricities of the ~ 700 merger binaries formed in run A300. The vertical dashed 
line is the hard-soft boundary for field stars of mass 1.0 Mq with one-dimensional velocity dispersion lOkm/s. The solid curves represent 
constant angular momenta J/ Jo = 0.2, 0.5, 1.0, and 2.0, where Jo is the total angular momentum of the system such that the pericentre 
of the initial hyperbolic orbit is 1.0 AU. 




a' (AU) 



Figure 13. Distribution of the scmimajor axes and eccentricities of the ~ 6000 merger binaries formed in run B005. The vertical dashed 
line is the hard-soft boundary for field stars of mass I.OA/q with one-dimensional velocity dispersion lOkm/s. 
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Figure 14. Histograms of semimajor axes of the merger binaries formed in runs A300 and BOOS, relative to the initial binary semimajor 
axis. 
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Figure 15. Histograms of eccentricities of the merger binaries formed in runs A300 and BOOS. The dotted lines represent properly 
normalized thermal distributions. 
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Figure 16. Cumulative distribution of pcriccntrc distance for tlic merger binaries formed in case A. The daslied lines, from left to right, 
correspond to a = 0.05, 0.1, 0.2, 0.5, 1.0 and 3.0^4(7 respectively. Each curve is equivalent to /t(^cp), the fraction of triple mergers as 
a function of the effective expanded radius of the first collision product, for a given a. 
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Figure 17. i?o.95) Ro.9i ^0.5i ^o.i and i?o.05 as a function of a for case A, where Rf is the value of i?cp at which = /. 
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Figure 18. Dependence of Rj/a on / = /y, with / < 0.5, for three different mass combinations. Solid triangles, open squares, and open 
triangles correspond to cases A, B, and C, respectively. 

Table 5. Fits for C in eg. llUI for the different mass combinations considered. 



Case moo (A/q) moi (Mq) mi (Mq) C 



A 


1.0 


1.0 


1.0 


1.61±0.01 


B 


1.0 


0.5 


1.0 


1.99±0.01 


C 


0.5 


0.5 


1.0 


2.40±0.02 


D 


1.0 


1.0 


0.5 


1.64±0.02 


E 


1.0 


1.0 


1.2 


1.78±0.02 



Table 6. Parameters of the binary-binary scattering experiments, including the mass of each star, rriij, the semimajor axis of each 
binary, a^, and the normalized cross sections for strong interactions and at least one collision to occur. 



run 


moo {Mq) 


moi (Mq) 


mio (A/q) 


mil (Mq) 


ao (AU) 


ai (AU) 


(^strong / ^oo \ 


'^"coll>l f Voo\^ 


7r(ao + ai)2 \ Vc ) 


7r(ao + ai)2 \ Vc I 


I 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 


0.62 ±0.13 


0.12 ±0.07 


II 


1.0 


1.0 


1.0 


1.0 


0.1 


0.1 


0.09 ±0.03 


0.04 ±0.02 


III 


1.0 


1.0 


1.0 


1.0 


1.0 


0.1 


0.090 ± 0.003 


0.0020 ± 0.0005 


IV 


1.0 


0.5 


1.0 


0.5 


0.1 


0.1 


0.12 ±0.03 


0.0099 ± 0.0004 
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Figure 21. Orbital parameters of four kinds of binaries formed in run I, with /oxp = 5. 
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Figure 22. Orbital parameters of four kinds of binaries formed in run II, with /c: 
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Table 7. Linear fits for the branching ratio of T/Q(> 2Mq) as a function of /oxp, where /T/Q{>2Af0) = ^/cxp + B. Also shown in the 
last column is the normalized cross section for the formation of triple-star /quadruple-star mergers with masses > 2Mq for /exp = 5. 



run 


A 


B 


'^T/Q(>2Mq)(/cxp - 5) /VocY 
7r(ao + ai)^ V f c / 


I 


0.007 


0.007 


0.011 ±0.002 


II 


0.0485 


0.14 


0.035 ±0.003 


III 


0.0182 


0.024 


0.00041 ± 0.00005 


IV 


0.0249 


0.0616 


0.034 ± 0.006 
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